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ABSTRACT 


vThe  convective  dispersion  of  a  solute  in  steady  flow  through  a  tube  is 
analyzed,  and  the  concentration  profile  for  any  Pec let  number  is  obtained  as  a 
convolution  of  the  profile  for  infinite  Peclet  number.  Close  approximations 
are  obtained  for  the  concentration  profile  and  its  axial  moments,  by  use  of 
orthogonal  collocation  in  the  radial  direction.  The  moments  thus  obtained 
converge  rapidly,  and  the  concentration  profile  less  rapidly,  toward  exactness 
as  the  number  of  collocation  po-  nts  is  increased.  A  two-point  radial  grid 
gives  results  of  practical  accuracy;  analytical  solutions  are  obtained  at  this 
level  of  approximation. 
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SIGNIFICANCE  AND  EXPLANATION 


Convective  dispersion  plays  an  important  role  in  many  processes  of 
chemical  reaction  and  separation.  Such  systems  are  commonly  analyzed  by  use 
of  a  radially-averaged  diffusion  equation^  This  approach  has  been  popularized 
by  the  simple  results  of  Taylor  (1953)  and  Aris  (1956)  for  the  cross-sectional 
average  concentration ,  and  perhaps  also  by  the  belief  that  the  radial 
variations  of  concentration  are  too  small  to  require  detailed  consideration* 
Subsequent  investigators  (e.g.  Gill  and  Sankarasubramanian  (1970),  Chatwin 
(1977),  De  Gance  and  Johns  (1978))  have  extended  this  approach  to  shorter 
times  and  to  other  boundary  conditions* 

The  foregoing  approach  is  not  easy  to  apply  to  reaction  or  separation 
processes  if  the  fluid  properties  vary.  nor  if  the  kinetics  or  equilibria  are 
non-linear*  Therefore,  in  this  paper  we  consider  an  alternate  methods  we 
solve  the  full  diffusion  equation  by  orthogonal  collocation  in  the  radial 
direction*  Non-reactive  systems  are  emphasized  here;  reactive  ones  will  be 
studied  more  fully  in  the  sequel  to  this  p  _ 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  of  this  report.  f  { 
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NEW  DESCRIPTIONS  OP  DISPERSION  IN  FLOW  THROUGH  TUBES: 
CONVOLUTION  AND  COLLOCATION  METHODS 

Janes  C.  Wanq  and  Warren  E.  Stewart 

INTRODUCTION 

The  collocation  procedure  described  here  is  a  variant  of  that  qiven  by 
Villadsen  and  Stewart  (1967);  see  also  Finlayson  (1972)  and  villadsen  and 
Michelsen  (1978).  For  axially  symmetric  states  in  a  tubular  reactor,  we 
approximate  the  profile  of  each  unknown  state  variable  Y^  in  the  form 


2i 


i=0 


lm 


0  <  £  <  1 


(1) 


This  expansion  can  also  be  written  as  a  Laqranae  polynomial 


n+1 


Ym  =  l  W  (C)Y  (T,Z,e  ) 


0  <  i  <  1 


k=1 


(2) 


with  radial  basis  functions 


wk(C)  =  n  a2  -  C2)/  n  (^  -  c2) 

j^k  3  jj*k  3 


(3) 


Thus,  each  approximate  profile  Y  is  represented  by  Laqranqe  interpolation 

m 

in  terms  of  its  values  Ym(TrZfC^)  at  a  set  of  chosen  radial  nodes  We 

2  2 

choose  £  ,  •••£  as  the  zeros  of  the  shifted  Legendre  polynomial  P  (x) 
In  n 

(Abramowitz  and  Stegun,  1972)  and  choose  £  as  the  wall  location,  K  *  1. 

n+i 

Any  needed  derivatives  or  integrals  of  Y  can  now  be  represented  as 

m 

linear  combinations  of  the  nodal  values  Y  (£.  )?  this  is  known  as  the  method 

m  k 

of  ordinates  (Villadsen  and  Stewart,  1967)*  For  example,  the  dimensionless 
radial  derivatives  of  interest  at  the  radial  node  5^  are: 


dY 


n+1 


'  m  ’ 


(4) 


n+1 

J,  \*WT'Z’ 
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.  ,  ax  1  ,  ,  uif, 

TdF  ^  dT^kz,^  *  J,  rId5  ^  dr)Ui1Ym(T'z'5k) 


n+1 

A  Bu,V(,'z>  • 


In  this  way,  we  reduce  all  radial  derivatives  to  summations;  the  same  can  be 
done  for  differences  and  integrals  with  respect  to  £•  For  the  above  choice 
of  nodes,  weiqhted  inteqrals  of  Y  over  the  cross-section  can  be  obtained 

m 

with  hiqh  accuracy  by  use  of  Gauss1  quadrature  formula* 

The  boundary  conditions  at  the  wall  may  be  linear  or  non-linear.  If  they 

are  linear,  then  Y  . „  can  be  expressed  linearly  in  terms  of  the  interior 
m ,  nri 

values  Y  by  use  of  the  boundary  condition  and  Eq.  (4).  This  permits 

ItIK 

elimination  of  Y  .  from  Eq.  (5);  the  resultinq  new  weiqht  coefficients 
m ,  n+i 

will  be  denoted  by  B^.  In  non-linear  cases,  on  the  other  hand,  one  retains 
the  wall  values  as  working  variables  and  applies  the  boundary  conditions 
numerically  during  the  solution  process. 


CONVOLUTION  RELATION 

Consider  the  developed  isothermal  flow  of  a  binary  Newtonian  or  non- 
Newtonian  fluid  in  an  infinitely  long  circular  tube.  The  continuity  equation 
for  either  species  may  be  written  in  the  dimensionless  form 


i?  +  v(5)  (5  S? 


an  i  a  aa,  i  a^a 


with  the  initial  and  boundary  conditions 


-  K1 1 '8 


At  T  .  0  :  ft(T,Z,e)  «  G(Z,£) 


At  5  -  1  :  &(T,Z,£)  +  K"ft(T,Z,S)  -  0  for  T  >  0 


i 


. - - -■  .  r..  ,  ■-  *•  i W  -WHHl 


a 

At  £  -  0  :  -yg-  Q(T,Z,e)  *  0  for  T  >  o  •  (9) 

Here  K,,#  and  K"  are  the  first  and  second  dimensionless  numbers  of 
Damkdhler  (1936)  for  homogeneous  and  heteroqeneous  reactions.  The  Laplace- 
Pourier  transform  of  Eqs.  (6)  -  (9)  from  (t,Z,£)  into  (s,p,£)  is 

2  «  H  jfj 

(a  ♦  k* 1  *  -  ^~j)n  ♦  pv{£)8  *  |  *5^*  (£  *5|*)  +  g(p,£)  (10) 


At  C  »  1  :  S(sfp#C)  +  K"B(8,p,£)  «  0  (11) 

At  5  *  0  :  !g>  S(s,p,C)  -  0  .  (12) 

Let  be  the  solution  of  Eqs.  (6)  -  (9)  with  Pe  =  Its  transform 

2 

satisfies  eqs.  (10)  -  (12)  with  the  p  term  suppressed*  Use  of  the  s-shift 
theorem  and  Fourier  convolution  (Campbell  and  Foster  1967)  then  gives: 

2 

fr(T,p,S)  -  exp(^-~.)^y.(T,p,C)  (13) 

Pe* 


0(T,Z,C)  -  r  “£S— e*Pl- -rlr  (2-X)2]Q_(T,X,C)dX  .  (14) 

-»•  2  /irr  4  T 

Thus,  the  complete  solution  of  Bps*  (6)  -  (9)  is  obtained  by  Gaussian 
smoothing  of  the  infinite-Pe  solution. 


TWO-POINT  COLLOCATION;  NON-REACTIVE  CASE 

Let  Cm  be  the  solution  of  Eqs.  (6)  -  (9)  with  K"  »  K* 1  *  *  0,  Pe  «=  °°, 
and  the  following  initial  mass  fraction  profile: 

G(Z,S)  -  «(Z)0(C)  .  (15) 

Collocation  with  n  *  2,  and  elimination  of  the  wall  concentration  throuqh 
the  boundary  condition  (8),  yields  an  initial-value  problem  for  the  nodal 


-3- 


functions  C0#(T,Z#?1  )  =  C1a>(T,Z)  and  C^r,Z,Z2)  =  C^{  T,Z).  Ihe  boundary 
condition  (9),  and  the  symmetry  of  the  true  solution,  are  satisfied  throuqh 
the  use  of  only  even  powers  of  £  in  the  collocation  basis  functions  fsee  Eq. 
(3)]«  In  place  of  Eq.  (10)  we  then  obtain 


sC 


i00 


+  pV 


k*1 


BJ  Z 

lk 


♦  0, 


1,2 


(16) 


in  which  and  pi  are  the  nodal  values  of  V(£)  and  $(£)•  To  complete 

the  solution,  we  evaluate  the  coefficients  from  Eqs.  (5)  and  (8),  with 

the  Gaussian  nodes 


52 .1.  JL 
1  2  m 


C2  =1  +  -j- 

2  2  /n  j 


(17) 


and  invert  the  transforms.  The  solution  differs  from  zero  only  within  the 
region  )z)  <  JtJ,  and  is  given  there  by 

0iP|t+z|i1 (X)  802Iq(X) 


C1<#  =■  exp(-8t  ){- 


2X 


+  ivf%r+V(z-T)} 


0  a|T-z I i  (X)  8Q  i  (x)  -  , 

C2«  -  exp(-8T){ - — -  +  |v~v^-  +  92«(Z+T)} 


(18.1 ) 


(18.2) 


in  which 


z  * 

Z  -  (V,  +  V2)T/2 

(19) 

T  - 

(V1  “  V2)T/2 

(20) 

6  ** 

256(V1  -  V2)"2 

(21  ) 

x  « 

/g(T2  -  Z2)  . 

(22) 

1 


& 


Interpolating  the  solution  accordinq  to  Eq.  (2),  and  intearatinq  over  the  tube 
cross-section,  we  find  the  mean  value  to  he  (C^  *  Gaussian 

quadrature  of  the  inteqral  qives  the  same  result* 


CENTRAL  MOMENTS  OP  THE  TOO- POT NT  SOLUTION 

Let  C  be  the  solution  of  Eqs*  (6)  -  (9)  and  (15)  with  K"  *  K' 1  *  *  0. 
The  mean  mass  fraction  in  a  cross-section  is  defined,  for  these  axisymmetric 
problems,  by 

<C>  =2  f1  CCd£  .  (23) 

The  central  axial  moments  of  <C>  are 

V.  =  /"  (Z  -  <V>t)i<C>dZ  (24) 

^  #  «pOP 

and  have  the  moment-qeneratinq  function 

Vi  =  (-  |^)i(exp(<V>Tp)  <C (T*p»S)>) | p_Q  (25) 


in  which  C(T,p,£)  is  the  Fourier  transform  of  C(T,Z,£)* 

For  the  collocation  solution  with  n  *  2,  we  thus  obtain 

"o  ■  7  '<>,  *  V 

»,  -  «s  -  v«*,  -  v’-'lP-' 


;  <VVgl  (i .  i-MPi-tiii,} 

M2  0X„  2  8  1 4  64  * 

Pe 


(26) 

(27) 


(28) 


For  a  radially  uniform  pulse  of  unit  strenqth,  p^  *  1;  and  for  developed 

2 

laminar  flow  of  a  Newtonian  fluid,  Vi  *  1  -  Eqs.  (26)  -  (28)  then  give 


Mo  =  1 


"1  w 

M2  =  (  2  +  ?92 1 2T 

Pe 


1-exp(-16T) 

1536 


(29) 

(30) 

(31) 
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The  exact  solutions  derived  by  Aris  (1956),  Gill  and  Sankarasubramanian 


(1970),  and  Chatwin  (1977)  for  a  radially  uniform  pulse  confirm  Eqs*  (29)  and 
(30).  Their  expressions  for  can  be  reduced  (after  a  stern  correction  in 

Aris*  solution)  to  the  common  form 


(T2 +  ik)2T 

pe 


1 


1440  +  32  t  Xj  exP<-XjT) 


(32) 


in  which  X^  is  the  jth  zero  of  the  Bessel  function  J^tx).  The  terms 
proportional  to  T  constitute  the  prediction  of  the  standard  Fickian 
dispersion  model,  with  the  long-time  asymptotic  dispersion  coefficient  found 
by  Taylor  (1953)  and  completed  by  Aris  (1956).  For  brevity,  the  latter  model 
is  called  the  Tfcylor-Aris  model. 

Before  testing  these  results  numerically,  we  give  collocation  solutions 
of  higher  order  for  the  axial  moments  of  the  function  C(t,Z,£). 


COLLOCATION  SOLUTIONS  OF  HIGHER  ORDER 

The  moments  of  the  mass  fraction  profile  about  the  origin  are  defined  by 
v*.«  •  c.  ZiJl(T,Z,C)dZ  .  (33) 

The  generating  function  for  these  moments  is 


1  p*0 

The  moments  at  finite  Pe  can  be  obtained  from  those  at  infinite  Pe  by  use 
of  Egs.  (13)  and  (34). 

A  collocation  solution  for  m^(T,$)  at  Pe  *  °*  is  obtained  as  follows. 
Let  be  the  vector  of  mesh-point  values  of  m^(x,£): 

We  apply  the  Fourier  transformation  and  n-point  radial  orthogonal  collocation 
to  Egs.  (6)  -  (9)  with  K"  *  Kf  • 1  *  0  and  Pe  *  Then  by  use  of  Eg.  (34) 
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’<* V *:i .  *  ■■ 


we  obtain 


It  Ti  *  B’»i  +  1  »  Hi-i  *  (36) 

Egs.  (13),  (34),  and  (36)  can  then  be  used  to  obtain  the  moments  at  finite 

values  of  Pe.  Finally,  the  cross-sectional  average  <mA>  can  be  obtained  by 

Gaussian  quadrature  of  . 

Tables  1,  2,  and  3  show  the  collocation  solutions  and  the  exact  solutions 
for  PQ,  U1*  and  for  two  initial  solution  distributions.  The 

collocation  solutions  converge  rapidly  toward  the  exact  results  for  all  values 
of  t,  and  for  both  initial  distributions.  From  Table  1  we  see  that  the 

collocation  solution  for  n  ■  2  is  particularly  accurate  in  the  case  of  a 

uniform  initial  radial  distribution. 

TWO-DIMENSIONAL  INITIAL  DISTRIBUTIONS 

Any  axially  symmetric  initial  distribution  can  be  expressed  in  the  form 

G(Z,5)  *  l  I  b  *j(Z)Qk<5)  (37) 

j  lc  3 

jJc  )c 

_  _  _  (2)  and  0  (£)•  C  (t,z,£)  be  the 

solution  of  Eqs.  (6)  -  (9)  with  constant  K"  and  K"',  and  with  the  initial 
distribution  cNo, Z,C)  *  6(Z)Qk(C).  Then  the  solution  corresponding  to  the 
initial  condition  (37)  is  obtainable  by  superposition, 

0(t,Z,5)  -  l  l  b  /"„  Vs  (X)Ck(T,Z-x,5)dX  (38) 

j  k  3 

and  has  the  Fourier  transform 

ff(T,p,C)  =  l  l  b  ^(pJcNt.p,^)  .  (39) 

j  k 

Let  and  Mk  be  the  ith  moments  of  the  functions  and  Ck(T,z,C) 

respectively  with  respect  to  Z*  Then  from  Bars*  (34)  and  (39)  we  get  the 
moments  m^  in  the  form 


Table  1  Second  moment  1000(11 - — )  for  initial  condition  Q(£)  *  1 

Pe 


T 

Collocation 

* 

Exact 

Taylor-Aris 

n=2 

n=3 

n=4 

.01 

.00791 

.00782 

. 00780 

.00779 

.1042 

.05 

.1623 

.1575 

.1574 

.1574 

.5208 

.10 

.5221 

.5059 

.5059 

.5059 

1.042 

.15 

.9705 

.9441 

.9442 

.9442 

1.562 

.20 

1.459 

1.425 

1.425 

1.425 

2.083 

.25 

1.965 

1.927 

1.927 

1.927 

2.604 

.40 

3.517 

3.474 

3.474 

3.474 

4.167 

1.0 

9.766 

9.722 

9.722 

9.722 

10.42 

* 

Exact 

first 

values  are  calculated  from 
two  moments  agree  for  all 

Eq.  (32).  The 
three  methods. 

Table  2  First  moment  100  for  initial  condition  Q(£) 


T 

Collocation 

★ 

Exact 

n=2 

n=3 

n=4 

.01 

-.1540 

-.1516 

-.1511 

-.1511 

.05 

-.5736 

-.5540 

-.5541 

-.5541 

.  10 

-.8314 

-.8083 

-.8086 

-.8086 

.15 

-.9472 

-.9298 

-.9299 

-.9299 

.  20 

-.9992 

-.9880 

-.9880 

-.9880 

.25 

-1.0226 

-1.0160 

-1.0159 

-1.0159 

.40 

-1.0399 

-1.0388 

-1.0388 

-1.0388 

1.0 

-1.0417 

-1.0417 

-1.0417 

-1.0417 

* Exact  values  are  calculated  from  Chatwin 
(1977).  The  first  moment  is  independent  of 
Pe,  in  view  of  Eqs.  (13)  and  (34). 


Table  3  Second  moment  1000(li  -  ~~z]  for  initial  condition  Q(£) 

_ Pe  _ 

* 

t  Collocation  Exact 


n=2 

n*3 

n=4 

.01 

.00791 

.00747 

.00738 

.00737 

.05 

.1623 

.1390 

.1389 

.1389 

.10 

.5221 

.4434 

.4447 

.4447 

.15 

.9705 

.8415 

.8431 

.8431 

.20 

1.459 

1.294 

1.296 

1.296 

.25 

1.965 

1.778 

1.779 

1.779 

.40 

3.517 

3.  305 

3.305 

3.305 

1.0 

9.766 

9.549 

9.549 

9.549 

* Exact  values  are  calculated  from  Chatwin 
(1977). 
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(40) 


■  U '■  V  1  ^<,'S> 

j  k  J  r-0 

for  any  initial  distribution  of  the  form  in  Eq.  (37)*  This  result  shows  the 
influence  of  the  initial  conditions  on  the  moments  of  the  mass  fraction 


profile  at  any  later  time. 


PROFILES  FOR  A  PULSE  OF  FINITE  LENGTH 

As  a  further  test  of  the  collocation  solutions ,  we  calculate  the  mass 
fraction  profiles  for  an  initial  pulse  of  maximum  amplitude  unity  and  of 
finite  lenqth  L  equal  to  0.1  unit  of  Z.  The  coordinate 

U  *  (22  -  T)/(L  +  T)  (41) 

is  introduced  here  to  facilitate  finite  element  calculations  of  the 
solution.  The  following  initial  condition  is  used, 

G(2,C)  =  *(U)  (42) 

in  which  v>(U)  is  a  C1  cubic  splines 


With  this  initial  condition,  the  solution  (valid  for  Pe  «  •)  is 

permanently  zero  outside  the  reqion  -1  <  U  <  1 • 

Fiqures  1 ,  2,  and  3  show  the  cross-sectional  averaqe  mass  fraction  <<*» 
as  a  function  of  t  and  z  at  Pe  «  15,  40,  and  infinity.  Three  methods  of 
solution  are  shown:  orthogonal  collocation  on  finite  elements  of  U,  two- 
point  radial  collocation,  and  the  Taylor-Aris  model.  The  radial  collocation 
was  done  in  all  cases  as  described  above;  the  finite  elements  of  U  were  C1 
cubic  splines  with  length  Au  *  0.025.  Each  method  was  first  applied  with 
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Pe  *  •,  and  Eq.  (14)  was  then  used  to  qet  results  at  finite  Peclet  numbers 


two-point  collocation  profiles  (dashed  curves)  model  the  true  solutions 
very  well,  except  for  some  lack  of  smoothness  at  short  times  with  Pe  » 

The  Taylor-Aris  model  exaqqerates  the  width  of  the  tails  of  the  pulse  in  every 
case. 


-  Reference  solution  by  finite  eleaents 

-  Collocation  solution  with  n  »  2 

»*n  Solution  by-  Taylor-Aris  nodal 


mi  Solution  by  Taylor-Aris  modal 


NOTATION 


coefficients  for  gradient  operator  in  Eg.  (4) 
aim'  bjk  coefficients  of  basis  functions. 

coefficients  for  Laplacian  operator  in  Eq,  (5) 

B’^  coefficients  for  Laplacian  operator  with  boundary  point 
eliminated. 

B*  matrix  of  order  n  with  elements 
m  ik 

C  solution  of  Eqs.  (6)  -  (9)  with  constant  K"  and  K* 1 1  and 

with  initial  condition  6(z)0(€) 

C*  special  form  of  C  for  initial  condition  i(z)QJc(€) 

limit  of  C  as  Pe  ♦  m 

G(5,z)  initial  distribution.  Eg.  (7) 

In(x)  modified  Bessel  function  of  the  first  kind 

first  and  second  Damkdhler  numbers 

m^  ith  axial  moment  of  ft  defined  in  Eq.  (33) 

m^  vector  of  nodal  values  of  m^r  Eg.  (35) 

n  number  of  collocation  points  interior  to  £  m  1 

p  Fourier  transform  variable 

Pe  *  RVmax^^AB  9  Peclet  number 

Q(£)  radial  function  in  initial  condition 

0^  value  of  p(5)  at  ith  collocation  point 
)( 

P  (?)  kth  radial  basis  function  in  initial  condition 
s  Laplace  transform  variable 
U  transformed  Z  coordinate  in  Eg.  (41) 

V(C)  velocity  profile,  relative  to  centerline  value 
value  of  V(C)  at  ith  collocation  point 
y  diaqonal  matrix  with  elements  vu  -  vi 


Y  mth  state  variable  (e*q*,  temperature  or  species  mass  fraction) 
m 

Z  “  z^B^Vmax  dimensionless  axial  coordinate 

<  >  averaqe  over  flow  cross-section 
(*)  binomial  coefficient 

Greek  Letters 

ith  central  moment  of  C  defined  in  Bq*  (24) 

£  *  r/R,  dimensionless  radial  coordinate 

t  *  dimensionless  time 

^ (U)  axial  function  in  initial  condition  of  Bq*  (43) 

*P^(z)  jth  axial  basis  function  in  initial  condition  of  Bq*  (37) 
a)  solute  mass  fraction 

fl  solution  for  u)  with  qeneral  initial  condition  G(z,£) 
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